In [1]:
import os
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import BeamDynamics as bd
import SimulationData as sd
import RFTrackTools as rfttools
import copy
import json
/home/tia/Repos/GIT_PSIPositronProduction/BdPyTools/BeamDynamics.py:11: UserWarning: ROOT module not available.
  warnings.warn('ROOT module not available.')
In [2]:
from importlib import reload
reload(bd)
reload(rfttools)
Out[2]:
<module 'RFTrackTools' from '/home/tia/Repos/GIT_PSIPositronProduction/RunningSimulations/RFTrack/RFTrackTools.py'>
In [3]:
# %matplotlib inline
# %matplotlib notebook
%matplotlib widget
plt.rcParams['figure.figsize'] = [9.6, 6.4]
defaultColorCycle = plt.rcParams["axes.prop_cycle"].by_key()['color']
# plotFont = {
#     'family' : 'sans-serif',
#     'weight' : 'normal',
#     'size'   : 12
# }
# matplotlib.rc('font', **plotFont)
# plt.rc('legend', fontsize=10)

Transverse Emittance Analisys¶

Distribution at Target Exit¶

In [4]:
DISTR_REL_PATH = '../../Data/Geant4/Pcubed/P3dist_baseline_p.ini'
targetDistr, _ = bd.convert_astra_to_standard_df(DISTR_REL_PATH, delim_whitespace=False, filterSpecsSelector=None, zProjection=None)
targetDistr.describe()
Out[4]:
Ekin Q betaRel gammaRel pdgId px py pz t x xp y yp z
count 137614.000000 1.376140e+05 137614.000000 137614.000000 137614.0 137614.000000 137614.000000 137614.000000 137614.000000 137614.000000 137614.000000 137614.000000 137614.000000 137614.0
mean 49.268378 1.000000e-09 0.995805 97.415810 -11.0 0.057750 0.002026 47.963867 0.092658 -0.007513 2.115502 0.002042 0.335502 0.0
std 118.740549 6.245236e-23 0.016748 232.369458 0.0 7.121371 7.055557 119.063066 0.006691 1.112848 448.462702 1.113585 447.292326 0.0
min 0.101143 1.000000e-09 0.550595 1.197931 -11.0 -87.036454 -131.671848 0.000602 0.000000 -13.017345 -1570.004841 -15.773611 -1570.700587 0.0
25% 7.736937 1.000000e-09 0.998079 16.140808 -11.0 -3.125669 -3.185303 6.016124 0.089999 -0.532107 -183.543793 -0.519669 -185.146241 0.0
50% 18.351462 1.000000e-09 0.999633 36.912915 -11.0 0.038935 0.003529 16.531034 0.092387 -0.004279 1.527649 0.002749 0.106735 0.0
75% 44.640179 1.000000e-09 0.999936 88.358651 -11.0 3.224203 3.175715 43.416008 0.094875 0.518339 188.208423 0.525416 185.503679 0.0
max 4972.376090 1.000000e-09 1.000000 9731.697275 -11.0 269.448714 82.801609 4972.880196 1.116313 16.176620 1570.451704 18.212887 1570.104422 0.0
In [5]:
plotSets = ['TransvPlane', 'TransvPsAngles', 'TransvPsMomenta', 'LongPsT']
plotDefs = bd.set_plot_defs_from_distrs([targetDistr], setNames=plotSets)
plotDefs[6]['lims2'] = [0., 100.]
_ = bd.plot_distr([targetDistr], plotDefs)
Figure
Figure
Figure
Figure
Figure
Figure
Figure

Overall Emittance¶

In [6]:
emitn = {}
emitGeom = {}
emitTraceSpace = {}
alphaTwiss = {}
betaTwiss = {}
gammaTwiss = {}
for planeName in ('x', 'y'):
    emitn[planeName] = bd.compute_emittance(targetDistr, planeName, norm='normalized')
    emitGeom[planeName] = bd.compute_emittance(targetDistr, planeName, norm='geometric')
    emitTraceSpace[planeName] = bd.compute_emittance(targetDistr, planeName, norm='tracespace')
    alphaTwiss[planeName], betaTwiss[planeName], gammaTwiss[planeName] = bd.compute_twiss(targetDistr, planeName)
    print(
        'emitn_{0:s} = {1:.1f} pi mm mrad, emitGeom_{0:s} = {2:.1f} pi mm mrad, emitTraceSpace_{0:s} = {3:.1f} pi mm mrad.'.format(
            planeName, emitn[planeName], emitGeom[planeName], emitTraceSpace[planeName]))
    print(
        'alphaTwiss_{0:s} = {1:.3f}, betaTwiss_{0:s} = {2:.3f} mm, gammaTwiss_{0:s} = {3:.3f} 1/mm.'.format(
            planeName, alphaTwiss[planeName], betaTwiss[planeName], gammaTwiss[planeName]))
Correcting offsets xAvg = -0.008 mm and pxAvg = 0.058 MeV/c.
Correcting offsets xAvg = -0.008 mm and pxAvg = 0.058 MeV/c.
Correcting offsets xAvg = -0.008 mm and xpAvg = 2.116 mrad.
Correcting offsets xAvg = -0.008 mm and xpAvg = 2.116 mrad.
Correcting offsets xAvg = -0.008 mm and xpAvg = 2.116 mrad.
emitn_x = 14981.7 pi mm mrad, emitGeom_x = 153.8 pi mm mrad, emitTraceSpace_x = 484.8 pi mm mrad.
alphaTwiss_x = -0.244, betaTwiss_x = 0.003 mm, gammaTwiss_x = 414.832 1/mm.
Correcting offsets yAvg = 0.002 mm and pyAvg = 0.002 MeV/c.
Correcting offsets yAvg = 0.002 mm and pyAvg = 0.002 MeV/c.
Correcting offsets yAvg = 0.002 mm and ypAvg = 0.336 mrad.
Correcting offsets yAvg = 0.002 mm and ypAvg = 0.336 mrad.
Correcting offsets yAvg = 0.002 mm and ypAvg = 0.336 mrad.
emitn_y = 14876.1 pi mm mrad, emitGeom_y = 152.7 pi mm mrad, emitTraceSpace_y = 484.3 pi mm mrad.
alphaTwiss_y = -0.240, betaTwiss_y = 0.003 mm, gammaTwiss_y = 413.096 1/mm.
In [7]:
emitRef = emitTraceSpace['x']

Transverse Emittance Vs. Bunch Portion¶

Values vs. 1D (?) Gaussian Distribution (To Be Revised)¶

In [8]:
gaussianPortions = {1: 0.6827, 2: 0.9545, 3: 0.9973, 4: 0.999937, 6: 0.99999998}
ellipseSpecs = {
    'x': {'alphaTwiss': alphaTwiss['x'], 'betaTwiss': betaTwiss['x']},
    'y': {'alphaTwiss': alphaTwiss['y'], 'betaTwiss': betaTwiss['y']}}
for planes in (['x'], ['x', 'y']):
    print('\nFiltering on {:d} planes:'.format(len(planes)))
    for Fa in (1., 2., 3., 4., 6.):
        selEllipseSpecs = {k: v for k, v in ellipseSpecs.items() if k in planes}
        distrWithinFaSigma, portionWithinFaSigma = bd.distr_within_ellipse(targetDistr, Fa**2.*emitRef, selEllipseSpecs)
        print(
            'Portion within {:.1f} sigma: {:.3f} (vs. {:.3f} for Gaussian).'.format(Fa, portionWithinFaSigma, gaussianPortions[Fa]))
Filtering on 1 planes:
Portion within 1.0 sigma: 0.596 (vs. 0.683 for Gaussian).
Portion within 2.0 sigma: 0.862 (vs. 0.955 for Gaussian).
Portion within 3.0 sigma: 0.958 (vs. 0.997 for Gaussian).
Portion within 4.0 sigma: 0.989 (vs. 1.000 for Gaussian).
Portion within 6.0 sigma: 0.998 (vs. 1.000 for Gaussian).

Filtering on 2 planes:
Portion within 1.0 sigma: 0.448 (vs. 0.683 for Gaussian).
Portion within 2.0 sigma: 0.793 (vs. 0.955 for Gaussian).
Portion within 3.0 sigma: 0.932 (vs. 0.997 for Gaussian).
Portion within 4.0 sigma: 0.980 (vs. 1.000 for Gaussian).
Portion within 6.0 sigma: 0.996 (vs. 1.000 for Gaussian).
In [9]:
refEllipseSpecs = {
    'x': {'alphaTwiss': alphaTwiss['x'], 'betaTwiss': betaTwiss['x']},
    'y': {'alphaTwiss': alphaTwiss['y'], 'betaTwiss': betaTwiss['y']}
}
FaList = np.arange(6., 0, -0.2)
distrWithinRef = []
bunchPortions = []
emitnPortions = []
for Fa in FaList:
    distrWithinRef.append(bd.distr_within_ellipse(targetDistr, Fa**2.*emitRef, refEllipseSpecs)[0])
    bunchPortions.append(distrWithinRef[-1].shape[0]/targetDistr.shape[0])
    emitnPortions.append(bd.compute_emittance(distrWithinRef[-1], 'x', verbose=False))

Emittance Vs. Bunch Portion¶

In [10]:
indsToAnnotate = [-5, -10, -15, -20, -25]
fig, ax = plt.subplots()
ax.plot(bunchPortions, emitnPortions, 'o-')
for ind in indsToAnnotate:
    ax.annotate(np.round(FaList[ind], 1), (bunchPortions[ind]-0.05, emitnPortions[ind]))
ax.set_xlim([0, 1])
ax.set_ylim([0, 14e3])
ax.set_xlabel('Bunch portion')
ax.set_ylabel('Norm. transv. emittance [mm mrad]')
ax.grid()
Figure

Export for Nico¶

In [11]:
# xy = np.column_stack([bunchPortions, emitnPortions])
# np.savetxt('emitnVsBunchPortion.dat', xy, header='{:24s} {:24s}'.format('bunchPortion', 'emitn [mm mrad]'))

Bunch Portion Vs. Cut Range¶

In [12]:
fig, ax = plt.subplots()
ax.plot(FaList, bunchPortions, 'o-')
ax.set_xlim([0, np.max(FaList)])
ax.set_ylim([0, np.max(bunchPortions)])
ax.set_xlabel('Fa of cut with Fa^2*rmsEmit ellipse')
ax.set_ylabel('Bunch portion')
ax.grid()
Figure

Visualization of Filtered Distributions¶

In [13]:
stepsToPlot = np.arange(5*2, len(distrWithinRef), 2)
FaListToPlot = [FaList[stepInd] for stepInd in stepsToPlot]
plotDefs[1]['lims1'] = [-15., 15.]
plotDefs[2]['lims1'] = plotDefs[1]['lims1']
plotDefs[1]['displayTable'] = False
ax = bd.plot_distr(
    [targetDistr, *[distrWithinRef[stepInd] for stepInd in stepsToPlot]], plotDefs,
    legendLabels=[
        'Full bunch',
        *['Within {:.1f} sigma'.format(Fa) for Fa in FaListToPlot]])
for ind, Fa in enumerate(FaListToPlot):
    bd.plot_ellipse(
        ax[1][0,0], Fa**2.*emitRef, semiAxisOrder=1, alphaTwiss=alphaTwiss['x'], betaTwiss=betaTwiss['x'], color='k')
    bd.plot_ellipse(
        ax[2][0,0], Fa**2.*emitRef, semiAxisOrder=1, alphaTwiss=alphaTwiss['y'], betaTwiss=betaTwiss['y'], color='k')
Figure
Figure
Figure
Figure
Figure
Figure
Figure

Transverse Emittance of Different Energy Windows¶

In [14]:
EkinWindowEdges = np.arange(0., 500., 2.)
particlesInWindow = []
emitnWindows = []
for wInd in range(len(EkinWindowEdges)-1):
    EkinWindow = {"Ekin": EkinWindowEdges[wInd:wInd+2]}
    particlesInWindow.append(bd.filter_distr(targetDistr, EkinWindow))
    emitnWindows.append(bd.compute_emittance(particlesInWindow[wInd], 'x', verbose=False))

Emittannce Vs. Energy¶

In [27]:
fig, ax = plt.subplots()
ax.plot(EkinWindowEdges[:-1], emitnWindows, 'o-')
ax.set_ylim([0, 15e3])
ax.set_xlabel('Ekin window bottom edge [MeV]')
ax.set_ylabel('Norm. emittance [pimmmrad]')
ax.grid()
Figure

Visualization of Filtered Distributions¶

In [16]:
windowIndsToPlot = np.arange(0, len(particlesInWindow), 1)
windowsToPlot = [EkinWindowEdges[wInd:wInd+2] for wInd in windowIndsToPlot]
ax = bd.plot_distr(
    [targetDistr, *[particlesInWindow[wInd] for wInd in windowIndsToPlot]], plotDefs,
    legendLabels=[
        'Full bunch',
        *['{:.1f} MeV < Ekin < {:.1f} MeV'.format(*EkinWindow) for EkinWindow in windowsToPlot]])
/home/tia/Repos/GIT_PSIPositronProduction/BdPyTools/BeamDynamics.py:1123: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
  fig, ax = plt.subplots(2, 2, figsize=(figWidth, figHeight))
Figure
Figure
Figure
Figure
Figure
Figure
Figure

Yield at DR Vs. Energy Window Size¶

Distribution at DR¶

In [17]:
DISTR_REL_PATH = {}
FILTER_SPECS = {}
REF_PARTICLE_PZ = {}
DISTR_REL_PATH['wp1'] = '../../Data/Astra/WP1atDR.dat'
FILTER_SPECS['wp1'] = {
    'z': [14480., 14530.],
    'pz': [0., 1620.]}
REF_PARTICLE_PZ['wp1'] = 1565.  # [MeV/c]
DISTR_REL_PATH['wp2'] = '../../Data/Astra/WP2atDR.dat'
FILTER_SPECS['wp2'] = {
    'z': [14470., 14520.],
    'pz': [0., 1620.]}
REF_PARTICLE_PZ['wp2'] = 1565.  # [MeV/c]
In [18]:
drDistr = {}
for wpName, relPath in DISTR_REL_PATH.items():
    drDistr[wpName] = pd.read_csv(relPath, names=bd.FILE_TYPES_SPECS['astra']['columnOrder'], header=None)
    # Minimal substitution of drDistr[wpName] = bd.extend_standard_df(drDistr[wpName], False)
    drDistr[wpName]['xp'] = 0.
    drDistr[wpName]['yp'] = 0.
    drDistr[wpName]['t'] = 0.
    drDistr[wpName]['Ekin'] = 0.
    drDistr[wpName]['pdgId'] = -11
    print('Working point {:s}:'.format(wpName))
    print(drDistr[wpName].describe())
Working point wp1:
             x        y             z       px       py            pz  \
count  30451.0  30451.0  30451.000000  30451.0  30451.0  30451.000000   
mean       0.0      0.0  14451.635141      0.0      0.0   1308.134190   
std        0.0      0.0    605.428368      0.0      0.0    422.126061   
min        0.0      0.0  -9991.133519      0.0      0.0      0.567927   
25%        0.0      0.0  14494.416345      0.0      0.0   1183.249370   
50%        0.0      0.0  14506.803113      0.0      0.0   1533.458841   
75%        0.0      0.0  14511.485757      0.0      0.0   1575.519662   
max        0.0      0.0  64988.153966      0.0      0.0   4883.870272   

         clock  macroCharge  particleIndex  statusFlag       xp       yp  \
count  30451.0      30451.0        30451.0     30451.0  30451.0  30451.0   
mean       0.0          0.0            0.0         0.0      0.0      0.0   
std        0.0          0.0            0.0         0.0      0.0      0.0   
min        0.0          0.0            0.0         0.0      0.0      0.0   
25%        0.0          0.0            0.0         0.0      0.0      0.0   
50%        0.0          0.0            0.0         0.0      0.0      0.0   
75%        0.0          0.0            0.0         0.0      0.0      0.0   
max        0.0          0.0            0.0         0.0      0.0      0.0   

             t     Ekin    pdgId  
count  30451.0  30451.0  30451.0  
mean       0.0      0.0    -11.0  
std        0.0      0.0      0.0  
min        0.0      0.0    -11.0  
25%        0.0      0.0    -11.0  
50%        0.0      0.0    -11.0  
75%        0.0      0.0    -11.0  
max        0.0      0.0    -11.0  
Working point wp2:
             x        y             z       px       py            pz  \
count  29279.0  29279.0  29279.000000  29279.0  29279.0  29279.000000   
mean       0.0      0.0  14459.170166      0.0      0.0   1421.584427   
std        0.0      0.0    366.967385      0.0      0.0    335.120112   
min        0.0      0.0 -11699.220033      0.0      0.0      1.290207   
25%        0.0      0.0  14492.560951      0.0      0.0   1501.087186   
50%        0.0      0.0  14496.303136      0.0      0.0   1557.978180   
75%        0.0      0.0  14497.844773      0.0      0.0   1574.104075   
max        0.0      0.0  21798.301070      0.0      0.0   3536.281605   

         clock  macroCharge  particleIndex  statusFlag       xp       yp  \
count  29279.0      29279.0        29279.0     29279.0  29279.0  29279.0   
mean       0.0          0.0            0.0         0.0      0.0      0.0   
std        0.0          0.0            0.0         0.0      0.0      0.0   
min        0.0          0.0            0.0         0.0      0.0      0.0   
25%        0.0          0.0            0.0         0.0      0.0      0.0   
50%        0.0          0.0            0.0         0.0      0.0      0.0   
75%        0.0          0.0            0.0         0.0      0.0      0.0   
max        0.0          0.0            0.0         0.0      0.0      0.0   

             t     Ekin    pdgId  
count  29279.0  29279.0  29279.0  
mean       0.0      0.0    -11.0  
std        0.0      0.0      0.0  
min        0.0      0.0    -11.0  
25%        0.0      0.0    -11.0  
50%        0.0      0.0    -11.0  
75%        0.0      0.0    -11.0  
max        0.0      0.0    -11.0  
In [19]:
for wpName, filterSpecs in zip(drDistr.keys(), FILTER_SPECS.values()):
    drDistr[wpName] = bd.filter_distr(drDistr[wpName], filterSpecs)
    print('Working point {:s}:'.format(wpName))
    print(drDistr[wpName].describe())
Working point wp1:
             x        y             z       px       py            pz  \
count  25337.0  25337.0  25337.000000  25337.0  25337.0  25337.000000   
mean       0.0      0.0  14506.027271      0.0      0.0   1414.660221   
std        0.0      0.0      7.661450      0.0      0.0    321.748731   
min        0.0      0.0  14484.702394      0.0      0.0      5.577844   
25%        0.0      0.0  14501.894891      0.0      0.0   1458.721448   
50%        0.0      0.0  14508.587808      0.0      0.0   1549.095748   
75%        0.0      0.0  14511.954643      0.0      0.0   1579.218798   
max        0.0      0.0  14527.074149      0.0      0.0   1619.837670   

         clock  macroCharge  particleIndex  statusFlag       xp       yp  \
count  25337.0      25337.0        25337.0     25337.0  25337.0  25337.0   
mean       0.0          0.0            0.0         0.0      0.0      0.0   
std        0.0          0.0            0.0         0.0      0.0      0.0   
min        0.0          0.0            0.0         0.0      0.0      0.0   
25%        0.0          0.0            0.0         0.0      0.0      0.0   
50%        0.0          0.0            0.0         0.0      0.0      0.0   
75%        0.0          0.0            0.0         0.0      0.0      0.0   
max        0.0          0.0            0.0         0.0      0.0      0.0   

             t     Ekin    pdgId  
count  25337.0  25337.0  25337.0  
mean       0.0      0.0    -11.0  
std        0.0      0.0      0.0  
min        0.0      0.0    -11.0  
25%        0.0      0.0    -11.0  
50%        0.0      0.0    -11.0  
75%        0.0      0.0    -11.0  
max        0.0      0.0    -11.0  
Working point wp2:
             x        y             z       px       py            pz  \
count  26491.0  26491.0  26491.000000  26491.0  26491.0  26491.000000   
mean       0.0      0.0  14494.991173      0.0      0.0   1489.906379   
std        0.0      0.0      5.332383      0.0      0.0    232.064203   
min        0.0      0.0  14472.672312      0.0      0.0      8.290856   
25%        0.0      0.0  14494.377212      0.0      0.0   1528.886562   
50%        0.0      0.0  14496.620400      0.0      0.0   1561.333039   
75%        0.0      0.0  14497.979419      0.0      0.0   1575.815748   
max        0.0      0.0  14519.997744      0.0      0.0   1598.811507   

         clock  macroCharge  particleIndex  statusFlag       xp       yp  \
count  26491.0      26491.0        26491.0     26491.0  26491.0  26491.0   
mean       0.0          0.0            0.0         0.0      0.0      0.0   
std        0.0          0.0            0.0         0.0      0.0      0.0   
min        0.0          0.0            0.0         0.0      0.0      0.0   
25%        0.0          0.0            0.0         0.0      0.0      0.0   
50%        0.0          0.0            0.0         0.0      0.0      0.0   
75%        0.0          0.0            0.0         0.0      0.0      0.0   
max        0.0          0.0            0.0         0.0      0.0      0.0   

             t     Ekin    pdgId  
count  26491.0  26491.0  26491.0  
mean       0.0      0.0    -11.0  
std        0.0      0.0      0.0  
min        0.0      0.0    -11.0  
25%        0.0      0.0    -11.0  
50%        0.0      0.0    -11.0  
75%        0.0      0.0    -11.0  
max        0.0      0.0    -11.0  
In [20]:
plotSets = ['z-pz']
plotDefs = bd.set_plot_defs_from_distrs(drDistr.values(), setNames=plotSets)
plotDefs[0]['lims1'] = [14488., 14516.]
plotDefs[0]['lims2'] = [1420., 1620.]
for distr in drDistr.values():
    _ = bd.plot_distr([distr], plotDefs)
Figure
Figure

DR Acceptance Vs. Energy Window Amplitude¶

In [21]:
def particles_in_window(distr, refParticlePz, pzWindowHalfAmplitudes, pzWindowTopEdgeScanRange, pzWindowAmplitudeStep):
    particlesInWindow = []
    for pzWindowAmplitude in pzWindowHalfAmplitudes:
        optYield = 0
        PZ_WINDOW_TOP_STEP = 1.  # [MeV/c]
        for pzWindowTop in np.arange(distr['pz'].max(), distr['pz'].max()-pzWindowTopEdgeScanRange*(1+PZ_WINDOW_TOP_STEP), -PZ_WINDOW_TOP_STEP):
            pzWindow = {
                "pz": [pzWindowTop-2.*pzWindowAmplitude*refParticlePz, pzWindowTop]
            }
            tmpParticlesInWindow = bd.filter_distr(distr, pzWindow)
            if tmpParticlesInWindow.shape[0] > optYield:
                optYield = tmpParticlesInWindow.shape[0]
                optParticlesInWindow = tmpParticlesInWindow
        particlesInWindow.append(optParticlesInWindow)
    return particlesInWindow
In [22]:
REF_WINDOW_HALF_AMPLITUDE = 0.038
MAX_WINDOW_HALF_AMPLITUDE = 0.06
PZ_WINDOW_AMPLITUDE_STEP = 0.001
PZ_WINDOW_TOP_EDGE_SCAN_RANGE = 200.  # [MeV/c]
pzWindowHalfAmplitudes = np.flip(np.arange(PZ_WINDOW_AMPLITUDE_STEP, MAX_WINDOW_HALF_AMPLITUDE+PZ_WINDOW_AMPLITUDE_STEP, PZ_WINDOW_AMPLITUDE_STEP))
particlesInWindow = {}
totParticlesInWindow = {}
for wpName in drDistr.keys():
    particlesInWindow[wpName] = particles_in_window(drDistr[wpName], REF_PARTICLE_PZ[wpName], pzWindowHalfAmplitudes, PZ_WINDOW_TOP_EDGE_SCAN_RANGE, PZ_WINDOW_AMPLITUDE_STEP)
    totParticlesInWindow[wpName] = [partInWin.shape[0] / drDistr[wpName].shape[0] for partInWin in particlesInWindow[wpName]]
    totParticlesInRefWindow = np.interp(REF_WINDOW_HALF_AMPLITUDE, np.flip(pzWindowHalfAmplitudes), np.flip(totParticlesInWindow[wpName]))
    totParticlesInWindow[wpName] /= totParticlesInRefWindow
In [26]:
fig, ax = plt.subplots()
for wpName in drDistr.keys():
    ax.plot(pzWindowHalfAmplitudes*100., totParticlesInWindow[wpName], '-')
    ax.plot(REF_WINDOW_HALF_AMPLITUDE*100., 1., 'ro')
ax.set_xlim([0., 5.])
ax.set_ylim([0., 1.2])
ax.set_xlabel('pz window amplitude [%]')
ax.set_ylabel('Particles in window')
ax.legend([wpName.upper() for wpName in drDistr.keys() for _ in range(2)])
ax.grid()
Figure

Export for Nico¶

In [24]:
# xy1y2 = np.column_stack([pzWindowHalfAmplitudes*100., *totParticlesInWindow.values()])
# np.savetxt('drAcceptanceVsPzWindowAmplitude.dat', xy1y2, header='{:24s} {:24s} {:24s}'.format('pzWindowHalfAmpl [%]', 'DRAcceptanceWP1', 'DRAcceptanceWP2'))

Visualization of Filtered Distributions¶

In [25]:
stepsToPlot = np.arange(0, len(particlesInWindow['wp1']), 10)
pzWindowHalfAmplitudesToPlot = [pzWindowHalfAmplitudes[stepInd] for stepInd in stepsToPlot]
for wpName in drDistr:
    ax = bd.plot_distr(
        [drDistr[wpName], *[particlesInWindow[wpName][stepInd] for stepInd in stepsToPlot]], plotDefs,
        legendLabels=[
            'Full bunch',
            *['Within +- {:.1f} %'.format(pzWindowHalfAmplitude*100.) for pzWindowHalfAmplitude in pzWindowHalfAmplitudesToPlot]])
Figure
Figure
In [ ]: